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Equilibrium solvation in quadrupolar solvents 

Anatoli A. Milischuk and Dmitry V. Matyusho\Q 

Department of Chemistry and Biochemistry, Arizona State University , PO Box 871604, Tempe, AZ 85287-1604 

(Dated: February 2, 2008) 

We present a microscopic theory of equilibrium solvation in solvents with zero dipole moment and 
non-zero quadrupole moment (quadrupolar solvents) . The theory is formulated in terms of autocor- 
relation functions of the quadrupolar polarization (structure factors). It can be therefore applied to 
an arbitrary dense quadrupolar solvent for which the structure factors are defined. We formulate 
a simple analytical perturbation treatment for the structure factors. The solute is described by 
coordinates, radii, and partial charges of constituent atoms. The theory is tested on Monte Carlo 
simulations of solvation in model quadrupolar solvents. It is also applied to the calculation of the 
activation barrier of electron transfer reactions in a cleft-shaped donor-acceptor complex dissolved in 
benzene with the structure factors of quadrupolar polarization obtained from Molecular Dynamics 
simulations. 
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I. INTRODUCTION 

The understanding of kinetics of chemical reactions in 
non-polar or, more generally, non-dipolar solvents poses 
the necessity to describe thermodynamioiiSi 3 ^^^ and 
dynamioLSiSiiSiii aspects of non-polar solvation.^? The 
overall solvation free energy in a non-dipolar solvent (zero 
permanent dipole) can be approximately separated into 
the contribution from the solute repulsive core (cavita- 
tion energy), the attractive dispersion solvation, and elec- 
trostatic contributions from induced dipoles and perma- 
nent multipoles starting from the quadrupole moment. 
The cavitation energy is often described by models con- 
sidering the free energy necessary to insert the solute 
hard repulsive core into the solvent pl&iii^iiSiil The dis- 
persion energy is a very significant part of solvation en- 
ergetics even in polar solvents^ It is normally modeled 
by site-site Lennard-Jones (LJ) interaction potentials, al- 
though LJ parameterization is often obscure for molecu- 
lar excited states. 

The positive cavitation free energy and negative con- 
tributions from dispersion and electrostatic interactions 
cancel each other in the overall solvation free energy 
which often constitutes only a small portion of each com- 
ponent. The cavitation and dispersion energies, however, 
cancel almost identically when the absorption and emis- 
sion solvatochromic shifts are subtracted to form the op- 
tical Stokes shift or when the solvent reorganization en- 
ergy of electron transfer (ET) is calculatedii^ The re- 
maining contribution from dispersion interactions, nor- 
mally associated with mechanical solvationji2iS£ does not 
typically exceed 100 — 200 cm -1 , which is small compared 
to usual values of the Stokes shift arising from dipo- 
lar and quadrupolar solvation^ The electrostatic com- 
ponent of solvation free energy in non-dipolar solvents 
arising from (partial) solute charges interacting with sol- 
vent quadrupoles is the focus of the present equilibrium 
solvation theory. 
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The electrostatic component of non-dipolar solvation, 
in particular quadrupolar solvation, has received little 
attention2iiMi£i2ii22i£i24 compared to the very extensive 
literature on solvation in dipolar solventsi£i2£i2& In con- 
trast to dipolar solvation, where dielectric measurements 
provide the basis for modeling thermodynamics and dy- 
namics of solvationjSSiSLSSiSLiffliiil there is no obvious 
method to extract the dynamic and thermodynamic re- 
sponse functions in non-dipolar solvents from existing 
data. This paper aims at bridging this gap by formulat- 
ing a microscopic solvation theory in terms of correlation 
functions of quadrupolar polarization of the pure solvent. 
In this paper, we limit the solvent correlation functions 
to the spatial domain thus gaining access to the equi- 
librium solvation thermodynamics. Once dynamic cor- 
relation functions become available, the theory can be 
extended to solvation dynamics. 

The experimental evidence on electrostatic quadrupo- 
lar solvation comes from the realm of steady-stateSi^SiS 3 . 
and time-resolved^ optical spectroscopy, resonance Ra- 
man spectroscopy^^ and from the kinetics of electron 
transfer (ET) reactions. 36 Spectroscopic chromophores 
and ET solutes are normally large molecules with com- 
plex molecular shape. Because of the relatively short 
range of interaction of molecular solvent quadrupoles 
with solute partial charges, it might be critical to in- 
clude the correct molecular shape of the solute into the 
formalism. We therefore sacrifice some accuracy of the 
modeling compared to potentially more accurate pertur- 
bation models for spherical solutes^ in order to introduce 
the molecular shape of the solute with atomic resolution 
combined with molecular charge distribution specified by 
atomic partial charges. It appears that this is the first 
theory of quadrupolar solvation approaching this level of 
detail in describing the solute. 

The paper starts with the formulation of the problem 
(Sec. |nj followed by the formal theory of solvation ther- 
modynamics (Sec. lIIl|) . The perturbation theory is given 
in terms of structure factors of fluctuating quadrupolar 
polarization in the pure solvent discussed in Sec. II VI The 
theory is tested on computer simulations of model ionic 
and dipolar solutes in Sec. Ivland is compared to exper- 
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FIG. 1: Charge transfer complex combining dymethoxyan- 
thracene unit for the donor (D) and a cyclobutene dicar- 
boxylate derivative for an acceptor (A) connected by a bridge 
(B)i— ■ The present theory is applied to the activated kinetics 
of charge separation, DB*A— >D + BA~. 

imental ET kinetics in Sec. IVII ET kinetics in a donor- 
bridge-acceptor cleft molecule referred to as complex 1 
(Fig. ^) has been extensively studied by Waldeck and 
Zimmt— The present theory is applied to the kinetic 
data in benzene used as a solvent. As a consistency test, 
the parameters used for ET rates in benzene are applied 
to ET in acetonitrile with the previously developed the- 
ory of polar solvationi2Li22i The overall good agreement 
between theory and experiment is reported. 

II. CONCEPTUAL FRAMEWORK 

Solvation in polar solvents is defined by the coupling of 
the field Eo (r) of the solute to the dipolar polarization of 
the solvent P(r). The interaction potential of the solute 
(subscript "0") with the solvent (subscript "s") is then 
a composite effect of this coupling integrated over the 
space occupied by the solvent 

vo.\P]=- f P(r) • E (r)dr. (1) 
Jn 

Here, the dipolar polarization is defined by the density 
of permanent dipoles rrij in the liquid 

P(r)=^m^(r-r i ), (2) 

3 

where the sum runs over N molecules of the solvent with 
center-of-mass coordinates Yj . In the linear response ap- 
proximation (LRA), the above interaction energy is sup- 
plemented by the Gaussian HamiltonianSi 

H P \P] = \l P(r) ■ X p(r, i-T 1 ■ P(r')drdr', (3) 



where the polarization response function xp(r,r') gen- 
erally depends on the shape of the field source. For cer- 
tain geometries, e.g. for a parallel plate capacitor, the 
dependence on geometry can be eliminated. The polar- 
ization induced in the solvent by the external electric 
field Eo(r) is then obtained by minimizing the functional 
v 0s [P] + H P [P] in P(r) to yield 

P(r)= f X p(r,r')-E (r')dr'. (4) 
Jn 

Molecular quadrupoles couple to an inhomogeneous 
electric field with the gradient VEo(r) 

v 0s [Q] = -(1/3) / Q(r) : VE (r)dr, (5) 
Jn 

where the quadrupolar polarization is 

Q(r)=^Q i «5(r-r J -) (6) 

3 

and Qj is the molecular quadrupole tensor. Similarly to 
Eq. ©, the Hamilt onian of the quadrupolar polarization 
in the pure solvent is given in the Gaussian form 

Hq[Q] = \\ Q(r) : Xg(ry)- 1 : Q(r')drdr'. (7) 

j n 

The minimization of t>os[Q] + Hq[Q] in terms of Q(r) 
then leads to 

Q(r) = i / X Q(r,r'):VE (r')rfr'. (8) 
Jn 

The definition of response functions xp(r, r') and 
XQ(r,r') incorporates the non-local solvent response af- 
fected by finite-range microscopic correlations of molec- 
ular dipoles and quadrupoles. These correlations are ne- 
glected in the continuum approximation 

XP, Q (r,r') = <5(r-r')xP,Q(r). (9) 

The continuum approximation, used in dielectric contin- 
uum models of dipolar— and non-dipola»22i2i solvation, 
significantly simplifies the calculation of the polar re- 
sponse function. In particular, the dipolar response func- 
tion can be related to the macroscopic dielectric prop- 
erties of the solvent through the macroscopic material 
Maxwell's equations. 

The Maxwell's equation for the overall electric field in 
the dielectric E reads^i 

V ■ (E + 4ttP) = 47177 + 471-V ■ (V ■ Q). (10) 

The material equations are closed by defining the dielec- 
tric displacement D = E + 47rP which is connected to 
the density of external charges p by the relation neglect- 
ing the quadrupolar density in the right hand part of Eq. 

V D = 4tt P . (11) 
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The dielectric displacement is in turn related to the over- 
all electric field through the static dielectric constant 
e s , D = e s E. In the case of a parallel plate capacitor, 
D = Eo and one gets 

XP = (e s - l)/47re s . (12) 

More complex geometries of the field source require solv- 
ing the Poisson equation with the boundary conditions 
defined by the dielectric constant and the shape of the 
dielectric. This procedure establishes the widely used 
dielectric continuum approximation for the dipolar re- 
sponse. 

Many of the advantages of the dielectric continuum 
approximation disappear when applied to quadrupolar 
(non-dipolar) solvation. The main problem is that the 
continuum quadrupolar susceptibility XQ does not come 
to the material Maxwell's equations and is, therefore, not 
directly related to any well-established experimental pro- 
tocol. Measuring quadrupolar susceptibility is still a non- 
trivial experimental problem^ In practical continuum 
calculations of the quadrupolar response, the quadrupo- 
lar susceptibility is obtained by fitting the calculated 
response to solvation free energies from spectroscopy^ 3 - 
One would alternatively want to have the quadrupolar 
susceptibility from properties of a pure quadrupolar sol- 
vent unaffected by all the complexities of treating spec- 
tral band-shapes of complex molecular solutes. In addi- 
tion, one needs to know the limits of applicability of the 
continuum approximation [Eq. 0] to relatively short- 
ranged quadrupolar interactions. The high directional- 
ity of quadrupolar forces make them unlike candidates 
for mean-field theories which are successfully applied to 
short-range but more isotropic dispersion forces^&^i All 
these considerations call for the necessity of microscopic 
theories of quadrupolar solvation. Once such a theory 
is formulated, approximate solutions, e.g. the continuum 
limit, can be obtained in a more controlled fashion. The 
formulation of such a theory is a goal of this paper. 

In order to develop a theory applicable to solutes of 
complex shape, we use a particular expression of the LRA 
to obtain the chemical potential of solvation [1q s . Within 
the LRA, [iQ S can be calculated as the second cumulant 
of the solute-solvent interaction potential 4 ^ 

- MOs = 08/2) (Sv 2 0s ) Q , (13) 

where (3 = l/k^T, is Boltzmann's constant and T is 
the temperature. The statistical average (. ..) is over 
the solvent configurations around a fictitious solute with 
the solute-solvent interaction energy vq s eliminated from 
the total interaction energy. The calculation of the aver- 
age (• ■ ■ )o over the solvent configurations in the presence 
of the repulsive core of the solute is a major challenge 
for the theory development. The solute expels the sol- 
vent from its volume and may significantly modify the 
statistics of solvent fluctuations either by altering the 
local density profile of the solvent or/and the statistics 
of orientational fluctuations of the solvent molecules. 39 



The expulsion effect propagates through the entire sol- 
vent changing substantially the solvent response func- 
tion in strongly dipolar liquids4£4I This strong effect of 
the solute on solvent statistics arises from the long range 
character of dipole-dipole interactions accounted for in 
dielectric models through the boundary conditions in the 
Poisson equation. 

The quadrupole-quadrupole interaction is much more 
short-range compared to the dipole-dipole interaction 
(oc 1/r vs cx 1/r 3 ). One may expect that statistics 
of orientational quadrupolar fluctuations is not signifi- 
cantly altered by the solute which only expels the solvent 
quadrupoles from its volume. In this approximation, one 
can switch from the statistical average (■ ■ ■ )o to the av- 
erage (...) over the statistical configurations of the pure 
solvent by replacing vq s {y) in Eq. (|T3*)) with vq s {t)6{t), 
where 8(r) is a step function equal to zero inside the so- 
lute and equal to one everywhere else. 

We will denote the component of the solvation chemical 
potential arising from fluctuations of quadrupolar polar- 
ization in homogeneous solvent as /j^ s (superscript "Q" 
stands for quadrupolar polarization). This component is 
expected to be the major contribution to the quadrupo- 
lar solvation chemical potential. However, in addition 
to expelling quadrupolar polarization from its volume, 
insertion of a solute into a liquid creates a nonuniform 
density profile represented by the solute-solvent pair cor- 
relation function ho s (r). This correlation function, highly 
specific to the solute shape and the thermodynamic state 
of the solvent, can be reliably calculated only for simple 
solute geometries. The component of /iq s associated with 
/io s (r) will be denoted as fj,Q S (superscript "D" refers to 
the local density profile) . The overall chemical potential 
of solvation is a sum of the long-range component due to 
quadrupolar orientational fluctuations ("Q" component) 
and the short-range component due to the local density 
profile ("D" component): 

Mo s =M? s +m£. (14) 

III. PERTURBATION THEORY OF SOLVATION 

The solute-solvent interaction potential in Eq. J5J) can 
be written either in Cartesian or spherical coordinates. 
We will use Greek indexes for the Cartesian projections 
and Latin indexes for the spherical projections^ In the 
Cartesian projections, the solute electric field gradient 
and the quadrupolar polarization are given as follows 

M a 

Mr) = -V a E p = -V a V p ]7T~~aT (15) 

a=l |r 

and 

3 

The sum (index a) in Eq. (|15fl runs over the M solute 
(subscript "0") atoms with coordinates and partial 
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charges qfi . In Eq. (|l(j|) , the summation is over TV solvent 
molecules with centers of mass at Yj relative to which the 
quadrupole tensor Q a /3j is defined as: 

K 

Q afid = (1/2) ^(r") 2 ^,/^ - M' (17) 



Here, the sum (index a) runs over the K atoms of the 
solvent molecule with coordinates r, + r?rf (f" = rj/rj) 
and partial charges q°-. 

In the spherical projections one gets 



I'Os 



J2 j ^2ro(r)Q2m,j( r ) rfr . 



(18) 



where the spherical projections of the quadrupolar polar- 
ization 



Q2m(r) = ^ Q2m,]S(r - Tj) 



(19) 



are given in terms of the spherical quadrupolar tensor 

K 



Q2m,j = E^°) 2y 2m(??) 



(20) 



In Eq. 1121 III . Y2m(?) is the spherical harmonic of the sec- 
ond order, —2 < m < 2. The distribution of molecular 
charge is illustrated in Fig. [5] on the example of the ben- 
zene molecule the first non-zero multipole of which is 
quadrupole. 




FIG. 2: The schematic representation of the benzene 
molecule: Vj is the coordinate of the center of mass, f" are 
unit vectors in the direction of partial charges q®. XYZ is 
the laboratory coordinate frame, the Z' axis of the laboratory 
system X'Y'Z' is along the unit wavevector k = k/fc. 



Substituting vq s from Eq. (|18|) into Eq. I|13|) and 

switching to the k-space we obtain 



8tt <W (2tt) 

m.n V / 



dk - 

j>2m (-k) (f>2n ( k ) S mn (k). 



. (21) 

Here, </>2m is the Fourier transform of 4>2m, <^2n is the 
Fourier transform of , and 



Q 2 = (2/3)Q : Q 



(22) 



is the rotational-invariant "effective axial quadrupole 
moment" i& The correlation function S mn (k) in Eq. I|21|) 
does not depend on the orientation of the wavevector be- 
cause of the rotational isotropy of the solvent: 



An 



(23) 



l ,3 



Since /ig s is invariant with respect to rotations of the 
coordinate system we first consider <\>2m and S mn in the 
coordinate system X'Y 1 Z 1 in which the wavevector k is 
directed along the Z'-axis (Fig.[3J). The functions of k in 
this coordinate system will be specified with the prime. 
The wavevector k introduces axial symmetry in the oth- 
erwise isotropic liquid of solvent molecules. The opera- 
tion of statistical average must therefore commute with 
the operation of rotation about the wavevector k 



exp 



S' mn exp (il z jj = exp (i(m - 71)7) S' m 



(24) 

where l z is the operator of rotation through the angle 7. 
The condition of invariance requires m = n leading to a 
simplified form of Eq. (|21|l 



q _ PpQ 2 



(Ik 



., (2tt, 

In Eq. I|25|) we represent 4> !m (k) as 

fi°(k) = 4>' (k), 

0' 1 (k) = i(&(k) + &(k)), 

0' 2 (k) = ~(& (k) + &(k)), 



4>' (k)S' m (k). (25) 



(26) 



where (m = — m) and 



y m (k) 



6 2m (-k)0* 2m (k) 

47T 



(27) 



Quadrupolar structure factors S' m (k) in Eq. (|2*5|) are 
defined as 



= S' Q (k), 
S n {k) = S[{k) + S[(k), 
S' 2 (k) = S' 2 (k) + S 2 (k), 



(28) 



■5 



where 




2 m 



AQ') 



2m, jt 



(29) 



Three quadrupolar structure factors S' m (k), grouped 
according to rotational symmetry, form a minimal set of 
correlation functions describing collective fluctuations of 
the quadrupolar polarization. Note that two structure 
factors representing uncoupled longitudinal and trans- 
verse dipolar polarization are sufficient to describe orien- 
tational fluctuations in dipolar solvents i^ 7 *^ Numerical 
calculations are more convenient to carry out in Carte- 
sian coordinates in which the chemical potential of solva- 
tion and the structure factors can be re-written as (Ap- 
pendix EJ 



Pos 



PpQ 2 



E 



^0 m (k)S m (fc). 



(27T) 



(30) 



Here, (f> m (k) the structure factors S m (k) are expressed in 
the form of rotation-invariant tensor contractions in Eq. 
(|A6|) and Eq. (|A5|> respectively. 

The density component of the solvation chemical po- 
tential in Eq. I|14|) is calculated by perturbation ex- 
pansion with the solute-solvent pair correlation func- 
tion h 1 ^} (r) corresponding to the distribution of the sol- 
vent around the repulsive core of the solute unaffected 
by the solute-solvent interaction potential vq s (reference 
system)^ 



(Pp/18) I dr/4°;(r)(0(r):Q 



(31) 



IV. QUADRUPOLAR STRUCTURE FACTORS 

The fc-dependent quadrupolar structure factors de- 
termine the microscopic spatial correlations of the 
quadrupolar polarization in the homogeneous sol- 
vent. Dipolar structure factors of model dipolar 
fluids 37 51 - 52 and fluids defined by site-site interaction 
potentialo 6 ! 53 ! 54 ! 55 ! 56 ! 57 ! 58 have been rather extensively 
studied in the literature. On the other hand, calculations 
of the quadrupolar structure factors have never been at- 
tempted before. Note that both dipolar and quadrupolar 
structure factors are unavailable from experiment, and 
liquid state theories and computer experiment are the 
only source of this information. We present here the re- 
sults of Monte Carlo (MC) simulations of structure fac- 
tors of quadrupolar hard-sphere fluids (simulation details 
are given in Appendix [Bj . This is followed by an ap- 
proximate analytical theory aiming at a fast algorithm 
applicable to calculations of solvation thermodynamics. 
In addition, the quadrupolar structure factors for ben- 
zene used in modeling ET reactions in Sec. IVlll were ob- 
tained from Molecular Dynamics (MD) simulations of the 
12-site nonpolarizable, rigid force field (Fig. [2J with the 



quadrupole moment Q = 8.63 DxA (simulation details 
are given in Appendix Q . 




FIG. 3: (a) Quadrupolar structure factors for fluids of hard 
spheres with embedded point axial quadrupoles: p* = 0.8, 
(Q*) 2 = 0.1 (1), 0.3 (2), and 0.5 (3). (b) (S°{k) - 1)5/8 
(dash-dotted line), -(S 1 (k) - 2)15/32 (dashed line), and 
(S 2 - 2)15/8 (solid line) with S m (k) from MC simulations 
at (Q*) 2 = 0.3. 

The quadrupolar structure factors S m (k) of a fluid of 
linear quadrupoles can be represented in terms of scalar 
products of unit vectors along the principle axes of the 
molecule and the direction of the wavevector k = k/fc. 
Adopting the notation of Ref. IH^ . = (e, • &j) — (e, • 
k)(ej • k) and Tj = (fij • k), one gets 

^(3lf -l)(3T 2 -l)e lk ^V 

ij I 




^^-(l-TfXl-T 2 )] e* 



(32) 



The axial-quadrupole structure factors S m (k) can also be 
expressed in terms of projections of the solvent-solvent 
pair correlation function on rotational invariants as fol- 
lows 



S°(k) 
S\k) 
S 2 (k) 

where 



-P 



(h 220 (k) 



~ h 222 {k) _^ (k) 

16 - 



= 1 
= 2 
= 2 

h 22l (k) = I ji{kr)h 22l {r)r 2 dr. (34) 
Jo 



|p [2h 22 \k)-h 222 (k) + ^h 22 \k)), 
\p (2h 220 (k) + 2h 222 (k)-^h 224 (k) j . 



(33) 
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In Eq. h mnl (k) is the Hankel transform of the 

projection h mnL (r) of the solvent-solvent pair correla- 
tion function on the corresponding rotational invariant , 59 
ji{x) is the spherical Bessel function. 60 




FIG. 4: Quadrupolar structure factors for the 12-site (rigid, 
nonpolarizable) benzene;— T = 298 K (solid lines) and T = 
342 K (dashed lines). 

Equations (|32|1 were applied to calculate S m (k) from 
NVT MC simulations with varying quadrupole moment 
((Q*) 2 = (0.1,0.3,0.5) in Fig.Ek). It turns out that the 
three quadrupolar structure factors can approximately 
be brought to one master curve by proper rescaling (Fig. 
|3Jd). This observation indicates that the long-range pro- 
jection h 22A (k) is the main component of S rn (k) sug- 
gesting a simple perturbation approach to the calcula- 
tion of S m (k). Taking the liquid without the quadrupo- 
lar interactions as a reference we expand h 22A {k) in the 
quadrupole-quadrupole interaction potential truncating 
the expansion by the first order term. Equation l|33[l can 
then be re-written in terms of the two-particle (super- 
script "(2)") perturbation integral as follows 

S°(k) = l-l2y q I^(ka,p*), 

S\k) = 2 + 16y q I {2) (k<r,p*), (35) 

S 2 (k) = 2-Ay q I^(ka,p*), 

where 

2ir 

y q = —PpQ 2 I° 2 (36) 

is the reduced density of solvent quadrupoles 5 " and a is 
the hard-sphere diameter. Note that the next perturba- 
tion term will result in the three-particle perturbation 
integral which we do not consider here. An improvement 
of the present description can be sought in terms of a 
Pade-truncated^ perturbation expansion for S m (k). 
The perturbation integral in Eq. I|35|l 

/oo 
dxgi°J(x,p*)u(kax)/x 3 (37) 

is defined in terms of the fourth order spherical Bessel 
function, ji{x), and the solvent-solvent radial distribu- 
tion function gis (x,p*), where x = r/a. For a fluid of 



TABLE I: Third-order polynomials a„(p*) in reduced density 
p* used in Eqs. and 



p 


ai,p 


a 2,p 


a 3,p 


a,4,p 





0.0095 


0.0901 


0.0971 


0.0620 


1 


-0.0776 


0.3580 


-0.4647 


0.5574 


2 


0.0800 


-0.5657 


1.0580 


-0.7940 


3 


0.0123 


0.1290 


-0.4370 


0.4810 



hard-sphere quadrupoles, the perturbation integral de- 
pends on ka and the reduced solvent density p* . It can 
be approximated by a series of spherical Bessel functions: 

4 

I^(ka,p*) = ^2a n (p*)j n (ka), (38) 

n=l 

where each a n (p*) is a third-order polynomial in reduced 
density: 

3 

a n (p*) = ^a n . p (p*) p . (39) 

p=0 

The fitting coefficients a n>p obtained by using the cor- 
rected Percus-Yevick radial distribution function for hard 
sphere fluids^ 2 , are listed in Table [I] The fit covers the 
range of reduced densities p* = 0.5 — 1.0. 




FIG. 5: Quadrupolar structure factors for fluids of hard- 
sphere axial quadrupoles from the perturbation expansion 
(Eqs. 1351 - 138H . solid lines) and from the MC simulations 
(dashed lines); p* = 0.8 and (Q*) 2 = 0.3. 

The analytical equations for the structure factors [Eqs. 
(|35H 1381) ] are compared to MC simulation results for the 
fluid of hard-sphere axial quadrupoles in Fig. [SJ The 
agreement is particularly good for (Q*) 2 < 0.3. Equa- 
tions (|35|I " H38|) are also used to define the structure fac- 
tors of the 12-site model of benzene (Fig. |2J approxi- 
mated by the axial quadrupole with the magnitude given 
by Eq. l|22|l . A reasonable agreement is obtained even by 
using the hard sphere radial distribution function instead 
of the actual pair distribution function of the reference 
potential of benzene (Fig.EJ. 
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k, A 



FIG. 6: Quadrupolar structure factors for the 12-site benzene 
from the perturbation expansion (Eqs. (I35H - 138L solid lines) 
and from the MD simulations (dashed lines); T=298 K, p* = 
0.982, and [Q*f = 0.45. 



V. MODEL SYSTEMS 

The theory is first tested on two model solutes serving 
as reference for many solvation studies: spherical ion and 
spherical dipole, both dissolved in an axial-quadrupole 
solvent. The formal theory is tested against the MC sim- 
ulations. Two types of simulations have been carried out. 
In the first set, we obtain the axial-quadrupole structure 
factors (Appendix^ required as input to the formal the- 
ory [Eq. I|30fl ]. In the second set, we directly calculate the 
LRA chemical potential of solvation [Eq. 1)13)) ] from fluc- 
tuations of the solute-solvent potential ( Appendix lB| . 



Ion 



For a point charge go inside a spherical cavity one has 
i> ' (k) = and the only nonzero component is given by 



^<V k ) = (47r g0 ) 2 jfjkRps 



(kR 0s ) 2 ' 



(40) 



where Ro s = (<tq + a) /2 is the distance of the closest ap- 
proach of the solvent (diameter a) to the solute (diameter 
(To). Equations (|3Tj|) and l|3*T|) then yield 

-PfjLos =y q (q* ) 2 [l Q (r 0s ,p*,y q )+I D (r 0s ,p*,y q )] , 

(41) 

where (q^) 2 = /3q 2 /a. The perturbation integrals in Eq. 
(|4*T]| are 



I Q (r 0s ,p*,y q ) 



I D {ros,P*,y q ) 



nr 0s Jo 



dxjf{xr 0s )S°(x), 



°°h$(x)dx 



(42) 



The perturbation integrals depend on the solvent reduced 
density p* — per 3 , the solvent quadrupolar density y q [Eq. 
({3"rJjl ]. and the solute-solvent size ratio ro s = Ros/&- They 
are tabulated as polynomials of p* and l/ro s in Sec. IV Cl 
Equation (|41|l is the microscopic perturbation solution 
for solvation of an ion. Below we will also consider two 



approximations to the complete solution: continuum ap- 
proximation and single-particle approximation. The con- 
tinuum limit for the solvation chemical potential can be 
obtained from the microscopic formulation by assuming 
that 4>°(k) changes much faster as a function of k than 
does the structure factor S°(k). When this is true, one 
can put S°(k) ~ S°(0) in Eq. g3J}. The density compo- 
nent disappears in the continuum limiljS with the final 
result 



= y q S°(0) 



Ml! 



(43) 



where the superscript "C" refers to the continuum limit. 
When correlations between the solvent dipoles are ne- 
glected by assuming h 221 = one gets (S°(0) = 1) 



* "|2 



-PPOs = Vq 



(go*) 



(44) 



where the superscript "S" refers to the single-particle ap- 
proximation. 



B. Dipole 

For a solute represented by a point dipole mo at the 
center of a spherical cavity of diameter oq, 2 (k) = 0, 
and two other components are 



0°(k) = (4tt) : 
^(k) = (Any 



,9(m -k) 2 j 2 (fci? ( 



()s) 



5 (kR 2 s ) 2 ' 
: 3(m 2 -(mo-k) 2 )j 2 (fcj?o s ) 
10 (kR 2 s ) 2 ■ 



(45) 



The total solvation chemical potential becomes 

-fipos = y q {ml) 2 [l Q (r 0s ,p*,y q ) + I D (r 0s , p* ,y q )] , 

(46) 

where (mj) 2 = /3m§/cr 3 and the perturbation integrals 
are 

2 f°° 

I Q {r Qs , P \y q ) = — / dxjl(xr 03 )(3S°(x)+S 1 (x)), 
7™o fl Jo 



I D (r as ,p*,y q ) = 5 



h i °J(x)dx 



(47) 



A polynomial approximation for dipolar I®' D (rQ S , p* , y q ) 
is given in Sec. IV CI 

The continuum and single-particle limits for dipolar 
solvation are 

-/W.=y<[3S°(0)+S l (0)]^ (48) 

Or, 



0s 



and 



-PPOs = Vq 



r 0s 



(49) 
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C. Perturbation integrals for (j,q 3 

Algebraic expressions for integrals I®(ro s , p* , y q ) in 
Eqs. I|42l) and H47|) can be obtained by using the pertur- 
bation expansion for the quadrupolar structure factors in 
Eq. (|35|) . The following integrals need to be calculated: 



dxjl{xr 0s )ji{xy) 



1536rt 



ir 2 as +y 2 ), 2r 0s >y (50) 



and 



I 12y- 



dxjKxr^j^xy) 



■(5r 2 s y 2 - 14r 4 s ), 2r 0s < y 



ny 



15360rg s 



(80rg a + 30r 2 s y 2 - 9y% 2r Qs > y 



?nrL 2r 0s <y, 



30y 5 



(51) 



where x — ha and y — r/a. 

The perturbation integral I®(ro s , p* 1 y q ) in Eq. is 
then 

I Q (ro s ,p*,y q ) = +y q Iion(ro s ,P*), (52) 



°'0s 



where 



Iion(r 0s ,p*) = 28r 2 , s I 2 ,8( r 0s, P*) - lO^^os, P*) 



1 



h ,-2(ro s ,p*) 



(53) 



Similarly, /^(r 0s , p* ,y q ) for the dipolar solute in Eq. 147|) 



is 



where 



I Q (ros,P*,y q ) = — +y q Idi P o\ c (ros,P*), (54) 



' OS 



28 



/dipole 0"Os,P*) = - —h, S^Os, P*) - 
r. 

-h,a{ros,P*) 



64r 0s 

In Eqs. (JEIi and (|55|) . 

h.n( r 0s,P*) = 

and 

h.nijOs.P*) 



,-2(^0s,P*) 
^2,2(r0s,P*)- 



2r 0a 



JO) 

dx 

2r 0s x 



(55) 



(56) 



(57) 



The numerical values of the perturbation integrals in Eqs. 
(|53|1 and i|55|) were fit to polynomials in p* and \/r$ s as 
follows 



li 



a n (p*) 



J I *\ \ " "n f 



n=4 



14 



' O.s 



^dipoio(? , o s ,p*) = 



a n (p*) 



n— 6,n^9 



' O.s 



where a n (p*) are third-order polynomials in p* 

3 

p=0 



(58) 



(59) 



(60) 



The fit covers p* ranging from 0.5 to 1.0 and ro s ranging 
from 0.8 to 2.4: the coefficients a n are listed in Table HTl 



D. Comparison to MC results 

MC simulations of solvation of spherical ions and 
dipoles have been carried out to test the formal the- 
ory. A solute is chosen as a hard sphere, (Jq/<j = 1.8, 
with charge, (q^) 2 = f3q 2 /<r = 15, or with point dipole, 
(toq) 2 = /3toq/(J 3 = 15, at the center. The initial config- 
uration in the cubic simulation box is created to accom- 
modate the solute at its center and iV solvent molecules 
with their size adjusted to keep p* — 0.8. The de- 
tails of reaction-field and Ewald sum corrections for the 
solute-solvent and solvent-solvent interaction potentials 
are given in Appendix iBl 

The results for solvation of the ionic solute are listed in 
Tab. IIIII Columns 2-4 give the variance of the the solute- 
solvent interaction potential from simulations. Simula- 
tions of ion solvation at N — 256, 500, and 864 show a 
noticeable size effect. Therefore, the infinite-dilution re- 
sult (column 5) was obtained by extrapolating the data 
at various N to N — > oo. Column 6 gives the theoreti- 
cal /io s . Its separation into the quadrupolar orientational 
component, Pq s , and the density component, p,Q S , is given 
in columns 7 and 8, respectively. The single-particle re- 
sponse p,Q S (column 9) turns out to be surprisingly close 
to pos from MC simulations. 

Results for solvation of the point hard-sphere dipole 
are shown in Tab. IIVI No dependence of the calculated 
quantities on the system size has been observed in this 
case for N > 500. Columns 2 and 3 give the average 
of the solute-solvent interaction energy and its variance. 

The near equality of these numbers supports the use of 

Q 



the LRAii* Column 5 gives the orientational part pfi 
from Eq. I|46|l . This component is about half of the over- 
all solvation chemical potential as is seen from the com- 
parison of column 5 to columns 2 and 3. The density 
component Pq s from Eq. (|4^|) adds to p® s to give the to- 
tal p,Q S in column 4 which is uniformly higher than the 
results of MC simulations. 
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TABLE II: Perturbation integrals from Eqs. I58I - 16UH . Numbers in columns indicate the coefficients b np in Eq. 1601 . 



V 


a,4,p 






a,7, p 


a S,P 


l9,p 








1 

2 
3 


-1/8 
1.618 
-5.240 
3.508 


0.061 
-15.631 
49.322 
-33.3 


1/64 
59.278 
-190.549 
130.4 


-1.341 
-118.0 
393.41 
-274.236 


4.284 
132.458 
-469.0 

334.8 


-5.611 
-83.681 
323.172 
-237.6 


3.39 
27.214 
-119.348 
90.924 


-0.777 
-3.447 
18.254 
-14.506 


P 


&6,p 






£JlO,p 


-^dipole 

asii,p 


Cll2,p 




1l4.p 



1 
2 
3 


-5/24 
-1.615 
4.530 
-2.614 


-1/8 
11.195 
-32.342 
18.942 


5/64 
-23.650 
70.892 
-42.622 


1/128 
87.281 
-268.563 
167.772 


3.781 
-179.677 
531.415 
-333.548 


-8.604 
176.884 
-492.183 
307.841 


6.739 
-88.443 
229.817 
-142.424 


-1.8 
17.9 
-43.5 
26.614 



TABLE III: Formal theory and MC results for ionic solute. 



(Q*f 


N = 256 


500 

f3 2 {Svl)/2° 


864 


CO 


-/W 




-Pfig h 




0.1 


0.23 


0.24 


0.24 


0.26 


0.30 


0.18 


0.12 


0.18 


0.2 


0.42 


0.87 


0.44 


0.48 


0.60 


0.35 


0.25 


0.37 


0.3 


0.59 


0.60 


0.61 


0.65 


0.86 


0.49 


0.37 


0.55 


0.4 


0.73 


0.76 


0.77 


0.84 


1.12 


0.63 


0.49 


0.73 


0.5 


0.87 


0.89 


0.91 


1.00 


1.36 


0.74 


0.62 


0.92 


0.6 


0.99 


1.02 


1.04 


1.15 


1.60 


0.88 


0.72 


1.10 



°MC simulations. 
6 Eq. dl) , 

c Single- particle solution, Eq. 1441 . 



In order to pin down the origin of the overestimated 
values, we compare our results with the Pade form of the 
perturbation expansion for a point dipole solute* (col- 
umn 7). This latter solution, which is in overall good 
agreement with simulations, can also be split into the 
orientational and density components (columns 8 and 9). 
The comparison of the orientational components of fiQ S in 
columns 5 and 8 and the density components in columns 
6 and 9 shows that it is the latter part that is overesti- 
mated in the calculation based on Eq. (|3 II) . Finally, the 
single-particle estimate, /j,q s (column 10), compares well 
with /j,q s obtained from the quadrupolar structure factors 
(cf. columns 5 and 10). The relatively high weight of the 
density component in dipole solvation makes the single- 
particle approximation less reliable than in the case of 
ion solvation. However, in overall, the single-particle for- 
mula works well for quadrupolar solvation. This obser- 
vation is consistent with the earlier notion by Ladanyi 
and Maroncell£ that the collective nature of the solva- 
tion response diminishes for higher multipoles and that 
non-dipolar solvation is dominated by single-particle re- 
sponse. We also note that the use of the quadrupo- 
lar structure factors from the analytical equations [Eqs. 
(|35f) - l|38[l ] gives results nearly identical (deviation < 1%) 
to the values obtained with simulated structure factors. 



VI. COMPARISON TO EXPERIMENT 

In this section we apply the present model to the calcu- 
lation of ET rates in the donor-acceptor complex shown 
in Fig. ^ The density component of fio s is hard to esti- 
mate for a molecule of such complex shape. On the other 
hand the charge in the charge-separated state D + BA~ is 
located close to the molecular surface. One might then 
expect that the solute-solvent interaction is locally of the 
ion-quadrupole type. Our calculations in Sec. II VI show 
that for this type of interaction potential the orienta- 
tional part constitutes about 75% of the overall solvation 
energy. The calculations below thus assume only the ori- 
entational component present in the solvent response 

Mo* - Vos- ( 61 ) 



A. Free energy gap and the reorganization energy 

Forward, fcf or (T), and backward, /cback(? 1 ), rates are 
available for complex 1 in benzene^^S^ whereas only for- 
ward rates have been measured in acetonitrile^i Rate 
measurements in benzene at different temperatures give 
access to the overall reaction free energy 



— /3A r G(T) = In (fc for (T)/£: back (T)) . (62) 
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TABLE IV: Formal theory and MC results for dipolar solute. 



(Q*) 2 -P(v 0s )/2 a (3 2 (Svl)/2 a -/W -P/ig b -I3^ s b -~P/i 0s c -f3)ig c -p^s c -Pvl 

0.1 0.59 0.52 0.52 0.26 0.26 0.51 0.27 0.24 0.28 

0.2 1.07 0.99 1.05 0.52 0.53 0.96 0.53 0.44 0.56 

0.3 1.48 1.41 1.56 0.77 0.79 1.38 0.78 0.60 0.84 

0.4 1.91 1.83 2.06 1.00 1.06 1.75 1.02 0.73 1.12 

0.5 2.24 2.19 2.56 1.24 1.32 2.10 1.25 0.84 1.40 

0.6 2.60 2.53 3.04 1.46 1.58 2.40 1.46 0.93 1.68 



°MC simulations. 
6 Eq. (gg). 

c Pade approximant from Ref. 0. The splitting into the orien- 
tational and density parts of the solvation energy is achieved by 
putting p* = in the perturbation integrals when the oricntational 
part is evaluated. 

"Single-particle solution, Eq. 1491 . 



TABLE V: Solvent parameters used in the calculations. 



Solvent 




Coo 


a, a A 


a, A 3 


m, D 


Q 


V 


, f> 


Benzene c 


1.00 


1.00 


5.28 


10.0 


0. 


8.63 


0.520 


544 


Benzene 


2.24 


2.24 


5.28 


10.4 


0. 


8.35 d 


0.515° 


544 


Acetonitrile 


35.0 


1.80 


4.14 


4.48 


3.9 


2.49 d 


0.424° 


100 



"From Ref. |63j 

6 In K, from Ref. lei 

Tarameters of the model benzene used in calculations. 
d In DxA, from Ref.0 



The free energy gap is a composite quantity including 
the vacuum energy gap A va cG, the quadrupolar free en- 
ergy A g G, the free energy of solvation by induced solvent 
dipoles Ai n dG, and the dispersion free energy Adi sp G: 



A r G — A vac G + A q G + A in( jG + Adi Sp G. 



(63) 



The difference in free energies is taken between the 
charge-separated final state ("f", D + BA~) and the ini- 
tial state ("i", D*BA) created by photoexcitation of the 
anthracene moiety of the donoriSiSSiSl 

The quadrupolar component of the energy gap is the 
difference of the final and initial values of the chemical 
potential of solvation 



A q G = 



AC (64) 
Using Eq. l(50jl . Eq. (|6*4"Jl can be re- written as follows 
(3 P Q 2 



A q G 



4tt 2 



/ dk(f™(k)-fr(k))S m (k), (65) 




k, A 



FIG. 7: Rotational projections fi"f,A from Eqs. H66II and I78H 
calculated for complex 1: m = (solid lines), m = 1 (dashed 
lines), and m — 2 (dash-dotted lines). Calculations are car- 
ried out for the charge distribution in the initial ("i") D*BA 
state (a), final ("f") D + BA~ state (b), and with the atomic 
charge distribution obtained as difference ( "A" ) of the atomic 
charges in the final and initial states (c). 



where 



= A; 2 <^>(k)) k . 



(66) 



The functions fj"f(k) are shown in Fig. 0("i" in (a) and 
"f" in (b)). 

The induction solvation is caused by the interaction 
of an induced solvent dipoles at point r with the elec- 



tric field of the solute. The free energy of induction sol- 
vation is obtained by integrating the real-space electro- 
static energy density Eij(r) 2 with the distribution func- 
tion go s (r) of the solvent molecules around the solute 



AindG 



(pa/2) / g 0s (r) [E f (r) 2 - E,(v) 2 ] dr. (67) 
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As in the case of quadrupolar solvation, we will neglect 
the solute solvent correlation function replacing <?os( r ) 
with a step function 9(r) which is equal to zero within 
the solute and is equal to one otherwise. By defining the 
Fourier transform of the electric field according to the 
relation 



E(k) = / E(r)(9(r)e ikr dr 



(68) 



one can rewrite Eq. I|67|) in the form of one-dimensional 
fc-integral 



/>oo 

AindG = -f^ / dkk 2 (S f (k) - Ei(k)). 



(69) 



Here, the density of electrostatic energy £(k) — £ L (k) + 
£ T (k) can be separated into its longitudinal and trans- 
verse components: 



^(fc) = (|(E(k).k)| 2 ) k 
£ T (fc)H|E(k)| 2 ) k -^(fc). 
The functions k 2 £ L ' T (k) are shown in Fig. [H] 



(70) 



6 

r 4 
2 


400 



•< 300 

x 

"-200 



100 


300 



= 200 
<] 



100 



A 1 1 

- 

/ ' 

" / ' V 


1 1 


1 1 

(a) - 


/ ' V 
/' V 

/' V 

II \\ 

— i V 

/' N>^-' 

/ 








1 1 


1 1 _ 

(b) - 




—J — 




rl 1 1 1 1 


1 1 


1 1 

(c) - 


_L / \ \ 
T i \ \ 
\i \\ 







0.5 



1.5 



2.5 



k, A 



FIG. 8: Longitudinal (solid lines) and transverse (dashed 
lines) projections of the electrostatic energy density [Eq. I7U> 1 
calculated for complex 1. (a) refers to the initial ("i") state, 
(b) refers to the final ( "f " ) state, (c) corresponds to the solute 
charge density obtained as difference ("A") in atomic charges 
in the final and initial states [Eq. I75H 1. 

The dispersion component Adi sp G is determined as the 
change in the total LJ solute-solvent interaction energy 



integrated over the solvent volume Q 



AdispG — p 



/ Au LJ (r)dr. 
Ju 



(71) 



Here 



UL.j(r) 



r-rn 



r-rn 



(72) 

£ as = v/eges and a as = (erg + cr)/2. The sum in Eq. 
(|72|l runs over all atoms in the solute. The atomic diam- 
eters, erg, and LJ energies, eg, are parametrized with the 
OPLS:§& The solvent parameters are listed in Table Ivl 

Calculations in acctonitrile were done without the 
quadrupolar component in the reaction free energy gap 
because of the very small reduced quadrupole moment 
of acetonitrile Q* relative to its reduced dipole m* (see 
Table F 1 ^ 



A r G — A vac G 



A P G 



AdispG. 



(73) 



The dipolar component, Ar,G, is calculated from the for- 
malism developed in Ref. |37] based on the integration of 
the longitudinal, x L (fc), and transverse, X T (&), compo- 
nents of the dipolar response function with the solute 
electric field: 



A D G = — 



W) 



, E 

P=L,T 



X p (k)£ p (k). (74) 



The reorganization energy of polar solvation X p is de- 
fined on the difference electric field Ea = E/ — E^ . The 
electrostatic energy density 



SaQh) = E A (k) 



(75) 



then separates into the longitudinal and transverse com- 
ponents resulting in corresponding components of the re- 
organization energy 

1 f dk 



Xr. 



(2tt) 3 E 

v ; P=L,T 



(76) 



In Eq. J7B), Xn T {k) 

are the nuclear response functions 
which, in contrast to x i:T (k), do not contain the effect 
of induced solvent dipolesi 37 ' 38 

The reorganization energy of quadrupolar solvation X q 
is given in terms of the gradient of E A (k). Following our 
formal derivation in terms of spherical coordinates, the 
final result can be written as 



PpQ 2 

4tt 2 



E 



dkf%(k)s m (k), 



where 



fZ(k) = k 2 (&(k)) l 



(77) 



(78) 



and A l (k) is obtained from the solute charge distribution 
build on difference atomic charges in the final and initial 
states (Fig.Et). 
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TABLE VI: Thermodynamics parameters (eV) of equilibrium solvation of complex 1 in benzene. 



T/K 


A q G a 


A q G 6 


A ind G 




Adis P G d 


A r G e 


A r G f 


A r G 9 


Ag 


X 6 


298 


-0.244 


-0.236 


-0.322 


0.196 


0.460 


-0.113 


-0.106 


-0.110 ft 


0.205 


0.198 


312 


-0.226 


-0.221 


-0.317 


0.193 


0.453 


-0.092 


-0.090 


-0.097 


0.189 


0.185 


326 


-0.218 


-0.208 


-0.312 


0.190 


0.446 


-0.082 


-0.085 


-0.082 


0.183 


0.174 


342 


-0.203 


-0.196 


-0.307 


0.186 


0.438 


-0.066 


-0.072 


-0.065 


0.170 


0.164 



°Eqs. JHSJ and g3 with S ' 1 ' 2 ^) from MD simulations. 
6 Eqs. JH3 and g3 with S°' 1 ' 2 (k) from Eqs. JHSJ and 1551 . 

c Calculated with Aa = 7.22 A 3 . 
Calculated with Aa = 17 A 3 . 

Calculated with Aa = 7.22 A 3 and A vac G = 0.258 eV. 
Calculated with Aa = 17 A 3 and A vac G = -0.68 X 10~ 3 eV. 
'Experiment from Eq. 1621 . 
Corresponds to T=296 K. 



B. ET rates 

The calculations according to the formalism outlined 
above require the combination of two input components, 
from the solute and from the solvent. The solute re- 
quires specifying coordinates, van der Waals radii (OPLS 
parameterization? 8 ), and charges of the solute atoms. 
The volume accessible to the solvent is limited by the 
solvent accessible surface (SAS) defined by complement- 
ing the solute atomic radii with the half of the solvent 
hard-sphere diameter, (erg + c)/2. Atomic coordinates 
and charges for complex 1 (Fig. [JJ are from Ref . The 
Fourier transforms of the solute field and solute field gra- 
dient are carried out on the 256 3 grid as described in 
Appendix iDl 

The calculation of the dispersion component AdispG 
requires knowledge of the alteration of the solute LJ en- 
ergy with electronic transition. Since this information 
is not available, we use the London dispersion potential 



to connect the change in the LJ energy to the change 
in the dipolar polarizability. The initial ET state is 
experimentally produced by photoexcitation of the an- 
thracene moiety of the donor. Anthracene polarizabili- 
ties in the ground (subscript "g" ) and excited (subscript 
"e") states areZi a g = 25 A 3 and a e = 42 A 3 . Accord- 
ing to the London equation e*, e ~ a ^/ e (" a " re f ers to 
an atomic site within the solute). The LJ energies on 
the anthracene atoms might hence be expected to scale 
as £° = (a e /a g ) 2 Eg leading to the scaling of the solute- 
solvent dispersion interaction potential as a e /a g . In the 
calculations below we will assume that LJ energies on 
only the anthracene moiety change with the transition. 
Since the anion acceptor state of complex 1 might in- 
volve some unknown polarizability change off-setting the 
polarizability change of the anthracene moiety, we con- 
sider Aa = a e — a g as a fitting parameter to reproduce 
the experimental A r G(T) [Eq. The other fitting 

parameter is A vac G. 



The second input to our calculation formalism comes 
from properties of the pure solvent. The polar, quadrupo- 
lar and dipolar, response comes into the theory in the 
form of three quadrupolar structure factors S m (k) and 
two dipolar structure factors S L - T (k). A parameteriza- 
tion scheme based on solvent dipole moment m, solvent 
diameter a, solvent polarizability a, and solvent number 
density p was developed previously in Refs. |3?] and |3j| 
Solvent parameters used in the calculation are listed in 
Table El 

The results of calculation of A r G(T) and X q (T) in ben- 
zene are listed in Tab. I VII Data in columns 2 and 10 in 
Table fVTI have been calculated by using S m (k) from MD 
simulations (12-site benzene). Equations l|35|) and i|38|) 
provide an alternative analytical route to S m (k) through 
the effective linear quadrupole moment with its magni- 
tude given by Eq. Results of using the analytical 
structure factors in Eqs. I|65l) and l|77|) agree well with 



corresponding values obtained by using the structure fac- 
tors from MD simulations (cf. columns 2, 3 and 10, 11 in 
Table Ejl. 

The reaction entropy AS* = —dA T G/dT is equal to 
— 8.55ke when Aa = 17 A of anthracene is used to 
calculate Adi sp G (column 6). The fit of experimental 
A r G(T) then yields A vac G = -0.68 x 10~ 3 eV. For these 
parameters, the dispersion and induction solvation en- 
tropies almost cancel each other, A^dS* = — 4.06ke 
(column 4) and AdispS" 3 = 5.80 kp (column 6), making 
the quadrupolar solvation entropy, A q 5° = — 10.29 ke 
(column 2), nearly equal to the experimental reaction 
entropy. The experimental value AS* = — 11.84 ke (col- 
umn 7) can be reproduced by downward scaling of Aa 
to 7.22 A 3 leading to A vac G = 0.258 eV. The dispersion 
solvation entropy then scales down to Adisp S° = 2.46 
(column 5). 

The parameters Aa and A vac G from the fit of experi- 
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TABLE VII: Rates (10 8 s" 1 ) of the ET in benzene. 



T/K 




^for 




h, , a 

^back 


. b 

">back 


r, c 
"-back 


298 


29.8 d 


28.54 


25.07 


0.48 d 


0.32 


0.36 


312 


27.5 


26.86 


24.64 


0.75 


0.79 


0.79 


326 


24.6 


26.02 


24.97 


1.31 


1.25 


1.11 


342 


21.5 


24.52 


24.78 


2.41 


2.45 


1.97 



"Experiment from Ref . l6fil 

'Theory with Ao=7.22 A 3 , A vac G = 0.258 eV, and \V\ = 8.0 
cm -1 . 

c Theory with A«=17 A 3 , A vac G = -0.68 x 10" 3 eV, and |V| = 
7.75 cm" 1 . 
d Measurement at T = 296 K. 



mental A r G(T) have been used to calculate ET rates in 
benzene (Tab. IVII|) according to the Jornter-Bixon for- 
mula (nonadiabatic ET): 



"^for/back 

(T) 



Tift 

a7 



1/2 oo 
n=0 



(79) 



exp 



P(\ S (T) ± A r G(T) + nhu v f 
4A S 



Here, S = X v /hv v is the Huang- Rhys factor and "+" and 
"— " refer to the forward ("for") and backward ("back") 
reaction rates, respectively. The vibrational reorgani- 
zation energy A„ = 0.39 eV and vibrational frequency 
hv v = 1410 cm -1 are from Ref. |HE The experimen- 
tal rates fcf or / back (T) are reproduced with the electronic 
coupling |V| = 8 cm -1 and A r G(T) and X q (T) from our 
calculations (columns 3 and 6 in Tab. IVIII see also Fig. 



moment presents an ideal choice. Table IVIHl lists the cal- 
culated thermodynamic parameters and rates of charge 
separation in complex 1 in acetonitrile in the experimen- 
tal temperature range between 255 and 335 K. 67 The 
calculations of the free energy gap and the reorganiza- 
tion energy are done by using the recently developed mi- 
croscopic theory of dipolar solvationi^S The dispersion 
part of the ET driving force in acetonitrile is calculated 
according to Eas. (|71|l and (|72[) with the solvent parame- 
ters listed in Table |V| 



TABLE VIII: Free energies (eV) and rate constants (10 s s 
for complex 1 in acetonitrile. 



T/K 


A P G 


AdispG 


A r G 


A p 




255 


-1.739 


0.101 


-1.381 


1.48 


2.91 


265 


-1.688 


0.099 


-1.331 


1.41 


3.28 


275 


-1.640 


0.098 


-1.285 


1.35 


3.66 


285 


-1.597 


0.096 


-1.243 


1.30 


4.03 


295 


-1.557 


0.095 


-1.204 


1.25 


4.39 


305 


-1.513 


0.094 


-1.161 


1.20 


4.64 


315 


-1.485 


0.092 


-1.134 


1.16 


5.07 


325 


-1.452 


0.091 


-1.103 


1.12 


5.40 


335 


-1.422 


0.090 


-1.074 


1.08 


5.70 





FIG. 9: Forward (circles and up triangles) and backward 
(squares and down triangles) ET rates in benzene for com- 
plex 1. Triangle symbols correspond to experimental datai— 
Circles and squares correspond to four temperatures at which 
5°* 1,a (fc) have been calculated from MD simulations. 



The parameters Aa and A vac G obtained from the fit of 
A r G(T) in benzene refer to the solute in the gas phase 
and do not depend on the solvent. The reliability of 
our fitting procedure can thus be further tested by us- 
ing these parameters to calculate ET rates in a different 
solvent. Since we also want a test independent of the 
present formulation for quadrupolar solvation, acetoni- 
trile with its small quadrupolc moment and large dipolc 



FIG. 10: The forward rates of charge separation in complex 
1 in acetonitrile. Open circles refer to experiment, the solid 
line refers to theory. 

Since the electronic coupling may depend on the 
solvent the experimental reaction ratesS 7 . were fit to 
Eq. (|79J) with electronic coupling V considered as the 
only fitting parameter fFig. I10|l. This procedure results 
in |V| = 2.21 cm -1 . This value is almost 4 times smaller 
than the corresponding electronic coupling in benzene. 
This trend parallels the one reported by Zimmt and 
WaldeckiSi V = 7.2 cm -1 in benzene vs V — 4.6 cm -1 
in acetonitrile. Qualitatively, this difference is attributed 
to the higher overlap of the donor and acceptor orbitals 
of complex 1 with the molecular orbitals of benzene re- 
siding in the clamp compared to the molecular orbitals 
of acetonitrile. 36 



VII. DISCUSSION 

Theories of dipolar solvation have been developed over 
the last decades to provide a hierarchy of approxima- 
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TABLE IX: Hierarchy of solutions for the solvation chemical potential, — /3/Xq s , in quadrupolar solvents. 



Approximation 


Solute 






Ion 


Dipole 


Single-particle 






Continuum 


(q5) 2 y q S° (0)/3r 3 3 


{m* ) 2 y q {3S°(Q) + S 1 (0))/5r 5 0s 


Perturbation 


(9o) 2 2/«(V 3r cL +y q Iion(r s,P*)) 


{ml ) 2 y q ( l/rfj s + y 9 Jdipoie (r 0s , p* ) ) 


Arbitrary solute 


Perturbation 


(P 2 P Q 2 /2) £ m / ^(k)S- m (fe)dk/(2^) 3 





tions ranging from the simple estimates of the Born and 
Onsager equations to more accurate Poisson-Boltzmann 
equation^ and up to microscopic liquid-state calcula- 
tions in terms of perturbation expansions'*^ integral 
equations^ or density functionatf^S theories. Noth- 
ing of comparable scale exists for quadrupolar solvation. 
Kim and co-workers2i*22i2ii2i have formulated a contin- 
uum theory of quadrupolar solvation of spherical solutes. 
The theory provides a fast calculation formalism, but 
suffers from the approximation of the spherical solute 
shape and empirical parameterization of the quadrupo- 
lar solvent susceptibility [Eq. J7J]. The integral equa- 
tion theories of the RISM familjiSiS avoid the problem 
of separate calculations of dipolar and quadrupolar sol- 
vation by considering the site-site interaction potentials 
automatically incorporating the infinite-order multipolc 
expansion. This advantage comes at the expense of the 
necessity to solve rather complex integral equations on 



one hand and to supply the charge-charge structure fac- 
tors of the solvent characterized by its site-site force field 
on the other. Force field parameterization is available 
only for a few solvents, and the calculations of the charge- 
charge structure factors requires extensive computations. 
Therefore, one would want a formalism for quadrupolar 
solvation that would incorporate the quadrupole moment 
of the solvent along with some experimentally available 
parameters as input in the microscopic solvation formal- 
ism. The formulation of such a theory is the result of 
this paper. It turns out that the quadrupolar structure 
factors of non-dipolar solvents can indeed be rather ac- 
curately calculated based on experimental input param- 
eters for the solvent: quadrupole moment, density, and 
the effective molecular diameter. Based on these struc- 
ture factors we have formulated a hierarchy of approxi- 
mations of increasing sophistication and accuracy listed 
in Table US 



As the crudest estimate we offer the single-particle ap- 
proximation which allows one to calculate the quadrupo- 
lar solvation energy from two solvent parameters, the 
quadrupolar density y q [Eq. and the solvent diam- 

eter a. The continuum approximation includes S m (0). 
However, since 5°(0) ~ 1 and 5 1,2 (0) ~ 2, the continuum 
limit is in fact very close to the single-particle result. The 
next approximation is based on perturbation integrals de- 
pending on the reduced solvent density p* = pa 3 and the 
reduced distance of the closest solute-solvent approach 
re = (To/2a + 0.5. The number of solvent parameters 
thus rises to three: y q , a, and p* . Finally, the calcu- 
lation of solvation free energy of a solute of arbitrary 
shape requires three angular projections of the auto cor- 
relation function of quadrupolar polarization (last line in 
Table ITX*|) . The problem which is still remain unresolved 
is the absence of an accurate algorithm for the density 
component of the solvation chemical potential [Eqs. 
(jT3)l ]. This component can be included for spherical so- 
lutes since the pair solute-solvent spherically symmetric 
distribution function can be then computed with suffi- 
cient accuracy. For instance, in the case of a spherical 
dipole, the Pade form of perturbation theories^ gives a 
good agreement with simulations. 
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FIG. 11: ln(fcf or ) vs 1/T for complex 1 is calculated in benzene 
with A vac G = 0.173 eV and in acetonitrile with A vac G = 
0.19 eV. In benzene X q and A r G: (a) are fixed at T = 320 
K; (b) assume full temperature dependence. In acetonitrile 
X q and A r G: (c) are fixed at T = 275 K; (d) assume full 
temperature dependence. 



Several effective medium approximations have been 
proposed in the literature to deal with the problem of the 
local density profiled They replace the solute-solvent in- 
teraction potential in the orientational component of the 
chemical potential of solvation (vo s 8(r) in Eq. with 
an effective solute-solvent coupling including the informa- 
tion about the solute-solvent density correlations (direct 
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solute-solvent correlation function in the surrogate&iSiS 
and density functional^ theories or vertex in the mode- 
coupling formulation 11 ). However, all such approxima- 
tions are heavily dependent on the detailed knowledge of 
the solute-solvent distribution function which is available 
only for solutes of simple shape. There is still no accept- 
able theory adequately addressing the calculation of the 
/Iq s components for solutes of complex shape. 

The present calculation for complex 1 highlights the 
importance of the induction and dispersion components 
of solvation in the equilibrium energy gap of ET reactions 
(Table IVI|) . The contribution of induction and disper- 
sion forces to the solvent reorganization energy are given, 
however, by a higher order of the perturbation theory, 
and they are normally much smaller than the correspond- 
ing components of the free energy gapi££ Quadrupolar 
solvation then becomes the most significant contribution 
to the solvent reorganization energy in non-dipolar sol- 
vents. In the crudest approximation, the quadrupolar 
reorganization energy is proportional to the quadrupolar 
density y q : 



22 



\ cx q 2 oP l3Q 2 /(R + a/2f 



for an ion and 



X q cx m^pPQ 2 / (R + a/2) 5 



(80) 



(81) 



for a dipolc. The present theory thus predicts a negative 
slope of X q vs T: 



(dXq/d^p^-XqiT- 1 +a P ), 



(82) 



where ap is the isobaric volume expansion coefficient. 
The slope of X q vs P is positive and is proportional to 
isobaric compressibility (3t- 



{dXq/dP), 



(3 T Xq. 



(83) 



The negative slope of the solvent reorganization energy 
has been obtained here for both quadrupolar ( Table IVT|) 
and dipolar (Table IVIIIfl solvents. This property thus 
seems to be a universal signature of ET reorganization 
in polar (dipolar and quadrupolar) solvents. This par- 
ticular dependence of the reorganization energy on tem- 
perature results in an observable effect in reactions with 
low activation barrier as was first noted in Ref. [TiJ The 
classical Marcus activation barrier for ET passes through 
a maximum when plotted vs 1/T (Arrhenius plot) at the 
point of activationless ET X q<p (T) + A r G(T) = when 
the negative slope of A 9jP vs T is applied. Figure ^2 il- 
lustrates this point by showing ln(fcf or ) vs 1/T calculated 
on complex 1 for which A vac G was adjusted to allow the 
point Ag iP (T) + A r G(T) = to fall in the experimental 
range of temperatures. The calculations with full tem- 
perature dependent X q . p (T) and A r G(T) are compared to 
the those under assumption X q>p = Const, A r G = Const. 
No maximum is seen in the latter case. We note also that 
intramolecular vibrations tend to mask the appearance of 
the maximum. 
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FIG. 12: Emission energy of ADMA 80 vs a q [Eq. (JH^] in 
quadrupolar solvents: ethylene (1), CO2 (2), benzene (3), 
tctrafluorobenzene (4), toluene (5), 1,4-dioxane (5). Emission 
ener gies , solvent parameters, and Ro = 4.32 A are taken from 
Ref. l8(l The dashed line is drawn as regression through the 
solvents excluding CO2. 



Although the full version of the theory is preferable for 
the analysis of spectroscopy and ET kinetics, the sim- 
ple single-particle approximation is useful in analyzing 
the qualitative trends with changing non-dipolar solvent. 
In particular, quadrupolar solvation results in the sol- 
vatochromic shift of optical transitions resulting in the 
change of chromophore's dipole state — > rxif 



hAi>j- 



-2a q (m f - m,) • nij 



where in the single-particle approximation, 



2tt(3 P Q 2 



5(Rq + a/2)' 



(84) 



(85) 



Figure ^| illustrates this trend based on the emis- 
sion energies of ADMA reported by Khajehpour and 
KauffmaniSSi Note that emission frequencies are affected 
by dispersion and induction solvation and the deviation 
of the supercritical CO2 from the linear trend may be 
traced to its lower polarizability resulting in lower dis- 
persion and induction stabilization energies. 

The present formulation, combined with the previously 
developed model for dipolar solvation^LS covers two ex- 
treme cases of polar solvation - purely dipolar and purely 
quadrupolar solvents. The formalism is based on the cor- 
relation functions of dipolar and quadrupolar polariza- 
tion as input and, therefore, can be applied to an arbi- 
trary dipolar or quadrupolar solvent. This approach has 
allowed us to address the problem of microscopic calcula- 
tion of solvation thermodynamics for solutes with atomic 
resolution of geometry and charge distribution. The ap- 
plication to real large donor-acceptor molecules shows 
encouraging results. Generalization of the theory will 
require considering dipolar-quadrupolar solvents, which 
will be a subject of future work. 
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APPENDIX A: DERIVATION OF EQ. \M 

Here we provide the derivation of the solute field 
gradient and quadrupolar structure factors in terms of 
rotationally-invariant contraction of Cartesian tensors. 
The Cartesian components of Q2m and 02m are 



Q20 — 
Q21 = 
Q22 = 



-in 



1/2 



Qzz, 



(Al) 



n . (Qxx Qyy ^Qxy) 1 

Z47T 



and 



5 

I An 



bxz + i4>yz) , 



(A2) 



loo 



2q \-^xx 4>yy ~\~ ^i^xy) 



Also, for negative to = -to, Q 2 .m = (-l) m Q2m: <h,m = 
{-l) m 4>* 2m . For the structure factors in the X'Y'Z' co- 
ordinates (k is collinear to Z' , Fig.0) one has 



S'°(k) 



iVQ 2 



^ ] Q zz,iQ zz.j^ 



Since k' = (0, 0, 1) in X'Y'Z' coordinates and Q' a g is a 
symmetric and traceless tensors one gets 



' ZZ,Z ^ ZZ ,j 



k'.Qj-k' k'.Q'.k' , 



^xz,iQxz,j Qyz,iQyz 



k' • Q ■ • Q'i ■ k' 



k'-Q^.k'j (k'.Q^.k'j, 

(Qxx,i ~ Qyy,i)(Qxx,j ~ Qyy,j) ^Qxy,iQxy,j 

= 2Q;:Q - 

4(k'.Q' r q;..k') + (k'.Q^.k') (k'-q^-k'). 

(A4) 

The tensor contractions in Eq. (|A4|) are invariant under 
rotations of the coordinate system. Therefore Eq. (|A3|) 
can be rewritten as 



NQ 2 

20 

3N& 



(k ■ Q, ■ k) (k • Q, ■ k 



i,3 



E 



k ■ ' ' k 



k-Qj-k) (k-Qj-k 

5 



3Ag^ 



2Qi : Qj - 4 ( k ■ Q, • Q 3 ;■ k 



S°(A) 
S\k) 

S 2 (k) 



(A5) 

where k = k/ |k|. Also the spherical components <fi m (k) 
entering Eq. (|30|l can be given in the rotation-invariant 
form as 

1 



E 



kQi-k kQ,k 



6°(k) = -|(k.0.k)| 2 , 
6 1 (k) = i((k-0*>-k)-|(k-0-k) 



<?(k) = ^ (20* : - 4(k- <j>* ■ <f> ■ k) + |(£ • <j> ■ k)f 

(A6) 

where e/> stands for <p a p (k). A numerical algorithm for the 
calculation of (j) a p (k) used in the calculations of complex 
1 is outlined in Appendix iDl 



S'\k) = 



20 



3AQ 2 



'z t {Qxz,iQxz,j Qyz,iQyz,j) e 



1,3 



& (k) — o N Q->(y^A(Qxx,i Qyy,i)(QxxJ Qyy,])~ 
^Qxy,iQxy,j^) e ' J j ■ 



(A3) 



APPENDIX B: DETAILS OF MC SIMULATIONS 

1. Quadrupolar hard-sphere fluids 

MC simulations have been carried out on a system of 
N = 500 hard sphere molecules with axial quadrupoles 
placed in a cubic simulation box. Periodic boundary con- 
ditions, minimum image convention, and the ratio of 0.5 
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between the distance of interaction cutoff and the size of 
the simulation box have been adopted^ The cutoff of 
the quadrupole-quadrupole interactions is corrected by 
the reaction field of continuum dielectric with the dielec- 
tric constant e' = oo. Simulation runs of the average 
length of 9 x 10 5 cycles with 4 x 10 5 cycles used to cal- 
culate S m (k) were performed at (Q*) 2 = (3pQ 2 /a 5 = 
(0.1,0.2,0.3,0.4,0.5,0.6) and p* = pa 3 = 0.8. 



2. Dipole solute in quadrupolar liquid 

We use the reaction field (RF) correction for the cut- 
off of long-range electrostatic interactions in MC simula- 
tions. The interaction energy of a dipolar solute with the 
quadrupolar solvent is the sum of the interaction energy 
u 0s{j) w hh the quadrupoles residing within the cutoff 
sphere and the RF correction terms uj^, and Uq s D ' sg1 



u 



where 



u 0s Q U) = r ja K f j0 • Qj ■ fjo) (m • fjo) - 2 (r j0 ■ Q 3 • m )] , 

(B2) 

and r,- = Tj - r , r j0 = r j0 /r j0 . In Eq. {BTJ, u^ Q;RF is 
the interaction energy of the solute with the polarization 
of the dielectric continuum (dielectric constant e') outside 
the cutoff sphere with the radius R c 



D_D;sclf 
l 0s 



(Bl) 



solvent quadrupole with the polarization of the dielec- 
tric continuum outside the cavity induced by the mth 
(m =/= j) solvent molecule is 



,QQ;RF 



(jm) 



2(Qj :Q m ) e'-l 



R 5 C Ze> 
The self energy of the jth quadrupole is 



Finally, uf (j) in Eq. (|B5|I is equal to u^(j) from Eq. 

For an axial quadrupole, Qva,j = 
(Q/2) (3e u> je a j — 5 ua ), where 5 va is the Kronecker 
delta symbol. The interaction potentials are given by 



nn.se\ir-\ _ (Qj : Qj) 

,DQ, 



(B7) 



(B8) 
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20 (ej 
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35 (©j • f m j 
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QQ;RF (jm) 

QQ;self^\ 



DQ:RF 
l 0s 



6(r j0 ■ Qj • m ) e' - 1 



3Q 2 e' - 1 
2i?s 3e' + 2 
3Q 2 e' - 1 
2i?s 3e' + 2 



3 (g^ • G m ) 1 



3e' 



(B3) 



(B9) 



DD;sett 
Os 



in Eq. (|B1|> is the energy of solute 



The term u 

self-polarization by the dielectric outside the cutoff: 

1 



l 0s 



2 / 



R 3 . 2e' + 1 ' 



(B4) 



The energy of the jth solvent molecule in the RF ge- 
ometry is 



u ss (j) 



m^j 



uf s Q (jm)+u 



QQ;BF 



(jm) 



-u 



QQ;sclf 



(J)+U% D (J) 

(B5) 



3. Ionic solute in quadrupolar liquid 

For ion solvation simulations we use the RF approxi- 
mation for the solvent-solvent interactions [Eq. i|B9(l ] and 
the Ewald sum (ES) approximation for the solute-solvent 
interactions. The solute energy is 



N 



UOs 



J2 u 0s Q U) + nZ CMi , (B10) 

3=1 



The quadrupole-quadrupole interaction energy is 



u ? S Q ti m ) = -t- (35(f r , 



X (f mj • Q 



m L mj j 



20 (f mj ' Qj ' Qm ' I* 



where M 0s C ' sclf is the self energy term and u® (j) is 
the interaction energy of the solute with the jth solvent 
molecule. The self energy is the interaction energy of the 
solute charge with its images in replicas of the simulation 
box 



•2(Qi : Q m )), 



m ' 1 mj ) 

(B6) 



CC:sclf 
l 0s 
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( S.EW 


7T 
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7T 2 


V 2 


h 6 




180 



where r mj - = r m - 



Y-mj ffrnj j &nd Qj • Qr; 



The interaction energy of the jth 



(Bll) 

where the finite-size correction £,ew = 1-41865 is accord- 
ing to Hunenberger and McCammonjS <7o is the hard- 
sphere solute diameter, and L is the side length of the 
cubic simulation box. 
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The interaction energy of the solute with the jth sol- 
vent molecules is (qo is the solute charge) 



f2{r j0 ) + -3-/1(^0) 



4gQ7T 

3L 3 



exp(-fc 2 /4K 2 ) cos(k • ryo)(k ■ Qj ■ k), 

|k|#0 

(B12) 



where 



fi(x) = erfc(Kx) + 



2kx exp(— x k ) 



/a(z) 



4k 3 exp(— x k ) 



(B13) 



and erfc(x) = 1 — erf(x), erf(x) is the error function. 60 
The solute-solvent interaction potential consists of two 
parts. The first is taken in the real r-space, while the 
second is taken in the inverted k-space. The convergence 
parameter k defines the number of replicas of the sim- 
ulation cell to be taken in the real space. Additionally, 
it defines the number of k-vectors taken in the inverted 
space sum. When this parameter is sufficiently large the 
calculation the real space sum is restricted to the original 
simulation box— 

Finally, the interaction energy of the jth molecule of 
the solvent is 



u s (j) 



'(j) 



u 



QQ\RF 



(j) + u? s Q -' sclt (j)+u% C (j), 



(B14) 



where the first three terms are given by Eq. (|B9(I . 
uff(j) = Mp s Q (j) [Eq. |B12|l ]. and the last term is given 
by Eq. (IBTlll . 



APPENDIX C: 



DETAILS OF MD SIMULATIONS 
OF BENZENE 



The MD simulations were carried out with the force 
field of 12-site benzene by Danten et al£L (TableEJ). The 
site-site interaction is given by the sum of the Lennard- 
Jones (LJ) and Coulomb interaction potentials: 



E 



ah 



As 



ab 



q q 



(Cl) 



where the LJ parameters are taken according to the 
Lorentz-Bertholet rules: e ab = \ / e a e b and a ah = (a a + 
<J b )/2. All simulations were done with the DL_POLY 
molecular dynamics packaged We run MD simulations 
in the temperature range from 298 K to 342 K with a 14 
K step. The timestep in each simulation is 5 fs. All MD 
simulation are 10 ns long. 

We used the Nose-Hoover— thermostat with the re- 
laxation parameter 0.5 fs. The proper choice of the sim- 
ulation timestep and the relaxation parameter ensures 



low drift in the total energy of about 0.3%. Cut-off 
distance of 12 A is used in calculations of LJ inter- 
actions. Ewald summation parameters used in calcu- 
lations are: (1) the convergence parameter k = 0.265 
A -1 ; (2) the maximum wavelength fc™^ 2 = 7 A" 1 . The 
simulation box is constructed to include 125 benzene 
molecules in a cube with the side length L = 26.46 A 
at T = 298 K to reproduce the experimental mass den- 
sity of benzene^ pu — 0.874 g/cm 3 . The side length 
is adjusted at each temperature to give the correct ex- 
perimental value for the isobaric temperature expansion 
coefficient^ a p = 1.14 x 10 _3 K _1 . 



TABLE X: Force field of benzene used in the MD simulation. 



Interaction site" 


a a /A 


e a x 10 2 /(kcal/mol) 


q a /e 


C 


3.473 


83.1 


-0.153 


H 


2.945 


12.5 


0.153 



«r cc = 1-393 A, 



fCH 



1.027 A. 



APPENDIX D: CALCULATION ALGORITHM 

FOR ^ a/3 (k) 

The gradient of the solute electric field (/) Q( g(r) is cal- 
culated on the 256 3 grid in r-space (Fig. I13fl as defined 
by Eq. (|15fl . In order to speed up the calculations, we 
split the calculation of the Fourier transform 4> a p (k) into 
two regions. In region f2 1; the Fourier transform is cal- 
culated numerically with a fast Fourier transform (FFT) 
algorithm— In region f^, the multipole expansion of the 
set of solute charges is used, and the Fourier transform 
is calculated analytically. The total <j) a j3 (k) is the sum of 
Fourier transforms from each region 



(k) = 4>al3; fii (k) + <i>af3- Q 2 (k) . 



(Dl) 



The analytical part of Fourier transform for an arbi- 
trary distribution of solute atomic charges gg with coor- 
dinates rg, r-g < R\ is given by 



M 



' 

Ri 



j„_l(fci?l) 



- (4*^ + f a 0a h)K-ii^) + KkpP^)) 



(D2) 



Here, R\ is the radius of the sphere enclosing volume f2i, 
j n (x) is the spherical Bessel function of the order 
P n {x) is the Legendre polynomial of order n, Pn(x) — 
dP n (x)/dx and P£(x) = d 2 P n (x) / dx 2 , x a = k • fg. 

The volume of integration for the numerical FFT is a 
cubic box of side lenght 



igrid — / X R D 



(D3) 



obtained as a multiple of the maximum extension of the 
solvent inaccessible cavity, given by the radius R ma ,x- The 
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L gnd 



FIG. 13: The schematic representation of the calculation ge- 
ometry. Qq is the solvent inaccessible cavity, fii (space be- 
tween the sphere of radius Ri and f2o) is the region of the 
numerical Fourier transform of <f) a 0(r) [Eq. 1151 1. ^2 is the 
region where Fourier transform of cf> a g(T) is taken analyti- 
cally [Eq. 1D2H . i?max is the radius of the sphere defining 
solute maximum extension and Ri — R m&K + a/2. The length 
igrid = 9 x ii m ax is the size of the r-space grid of dimension 
256 x 256 x 256; (m, n, k) is the position of a point on the 
grid. 



radius of the sphere of the solute maximum extension is 
defined as 

fl max = max{|rg - r gc | + a«/2} , (D4) 

where i runs over M solute atoms. The sphere is cen- 
tered at the geometric center of the solute molecule, 
r gc = (1/M)J2 rg and Cq is the diameter of atom a 
of the solute. The radius of the spherical region Sli is 
the obtained by adding the solvent radius a/2 to i? m ax: 
Ri = -R m ax + cr/2 (Fig. I13|) . The multiplication factor 
/ = 9 in Eq. i|D3|) is chosen as a trade-off between mini- 
mizing FFT errors from artificial periodicity of the lattice 
sum and the need for a sufficiently small increment in the 
k-space for the k-integration. 

Charges of complex 1 in the non-polar (initial) 
and charge-transfer (final) states are calculated with 
Gaussian'03 (UHF/6-31G(3p))S Partial atomic charges 
are obtained by fitting to the electrostatic potential given 
by the exact wavefunction in the CHELPG schemed 
Charges are calculated separately for the isolated donor, 
donor cation, acceptor, acceptor anion, and the neutral 
bridge. The charges of the two hydrogens that substi- 
tute the rest of the clamp in the model compounds are 
incorporated in the connected carbon atoms^i 
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